Loss of least-loaded class in asexual 
populations due to drift and epistasis 



Kavita Jain 

Theoretical Sciences Unit and Evolutionary and Organismal Biology Unit, 
Jawaharlal Nehru Centre for Advanced Scientific Research, 
Jakkur P.O., Bangalore 560064, India 

May 30, 2008 



1 



Running head: Loss of least-loaded class 

Keywords: asexual evolution, Muller's ratchet, epistasis 

Corresponding author: 
Kavita Jain, 

Theoretical Sciences Unit and Evolutionary and Organismal Biology Unit, 
Jawaharlal Nehru Centre for Advanced Scientific Research, 
Jakkur P.O., Bangalore 560064, India, 
j ain@ j ncasr . ac . in 

Abstract: We consider the dynamics of a non-recombining haploid popula- 
tion of finite size which accumulates deleterious mutations irreversibly. This 
ratchet like process occurs at a finite speed in the absence of epistasis, but 
it has been suggested that synergistic epistasis can halt the ratchet. Using 
a diffusion theory, we find explicit analytical expressions for the typical time 
between successive clicks of the ratchet for both non-epistatic and epistatic 
fitness functions. Our calculations show that the inter-click time is of a scal- 
ing form which in the absence of epistasis gives a speed that is determined by 
size of the least-loaded class and the selection coefficient. With synergistic 
interactions, the ratchet speed is found to approach zero rapidly for arbitrary 
epistasis. Our analytical results are in good agreement with the numerical 
simulations. 
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In an asexual population of size N, even the fittest individuals can be 
lost by stochastic fluctuations arising due to the finiteness of the popula- 
tion size. If the individual's genome is long enough that the back mu- 
tations can be ignored and recombination is absent, the minimum num- 
ber of deleterious m utations ( l east-l o aded class) in a finit e population in- 



creases irreversibly (IMullerI Il964j ; IFelsensteini [1974]). For this rea- 



son, this process has be en invoked as a potential cause for the evolution of 



sex and recombination (IJudson and Normarki . 1 19961 ; IHurst and Peck . 



19961 : IBarton and CharlesworthI. Il998l). deg e nerat ion of nonrecombin- 



ing pa rts like Y chr omosome (ICharlesworthi . Il978l ) and mitochondrial 
DNA ( LoeweI . 2006 ) of sexually reproducing organis ms and extinction of ob - 
ligately asexual populations by mutational meltdown ([Gabriel et a/.l . ll993l ). 

Due to the irreversible accumulation of deleterious mutations, the pro- 
cess described above acts like a ratchet each click of which corresponds to the 
loss of the least-loaded class. In the simplest model known as the Muller's 
ratchet, the ratchet clicks at a constant rate which depends on the population 
size N, mutation rate U and selection coefficient s. The ratchet speed is also 
known to depend on other biologically relevant fact ors such as recombina- 



tion rate (IBellI [19 88; C harlesworth et a i. 1993 ), epistatic inter actions 



( Kondrashov . 1 1994 IButcherI . Ii995 ~ Ischultz and Lynch. Il997), frac- 



tion and selection coefficient of favorable mutations (IWOODCOCK and HlGGS 



199 61; IBachtrog and GordoI . 120041 ) and spatial structure of the population 



( ICOMBADAO et all 120071 ). Although much numerical data for the ratchet 



speed is available, very few analytical results are known. 

As it is desirable to stop or at least slow d own the ratchet, several mecha- 

nisms with this objective have been proposed (IBellI . Il988l : I Wagner and Gabriel! . 



19901 ; ICharlesworth et all Il993l) . One such proposal is to inc l ude epistatic 



inter actions in the genome fitness (ICharlesworth et a/.l . ll993l ; lKoNDRASHOvl . 
19941 ). It has been observed experimentally that the gene loci do n ot al- 



ways contribute independentl y to the genome fitness (IWolf et all . l2000l : 



DE Visser and ElenaI . |2007| ) and the effect of two deleterious mutations 



can be better (antagonistic) or worse (synergistic) than were they to act 
independently. For Muller's ratchet operating under epistatic selection, it 
was noted using numerical simulations that "sufficiently st rong synergistic 
epist asis can effectively halt the action of Muller's ratchet" (IKondrashovI . 
1994J). However it was not shown how the ratchet speed approaches zero 



asymptotically and how much epistatic interaction is required to halt the 
ratchet. 
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In this article, besides the classic Muller's ratchet that assumes haploid 
asexual population evolving on a fitness landscape in which each gene locus 
contributes independently to the fitness of the genome, we also study Kon- 
drashov's model which considers fitness functions with epistatic interactions. 
We assume that an individual with k mutations has a fitness 



W(k) 



:i-s) 



(i) 



where s is the selection coefficient. For a = 1, the epistatic interactions are 
absent while a > 1 corresponds to synergistically epistatic fitness. Our main 
purpose is to obtain explicit analytical expressions for the typical time Tj 
elapsed between the (J — l)th and Jth click of the ratchet in these models. 
If the time spent between any two successive clicks is found to be constant 
then the ratchet turns with a finite speed 1/T, while it is said to be halted at 
large times if the inter-click time increases with the number of accumulated 
mutations. 

In the past, Muller's ratchet has been investigated using a diffusion ap- 
proximation which assumes that the popul ation n n of the least loaded class is 



large and applies to slowly c l icking ratchet (IStephan et all Il993l ; IGORDO and Charlesworth 



2000al ; IStephan and Kimi . 12003 ) . The opposite situation of small n n and 



high r atchet rate has been described by a momen t s method (IGabriel et al. 



1993; IGesslerI 1199 5; Hi ggs and Woodcock! . Il995l ; IPrugel-Bennett . 



19971 ; IRouzine et all 120031 ) . In this article, we adopt the first method that 
works in the parameter range for which no 3> 1 and the ratchet clicks are slow 
enough that the population can equilibrate between successive clicks. This re- 
quires the knowledge of the steady state properties of a n infinitely large popu- 
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, 1994) 



but have been studied numerically for epistatic case 



we solve the deterministic quasispecies equation in steady state for a > 1 and 
show that the population frequency of the class with minimum number J of 
mutations increases with J. These deterministic results are then used to find 
an expression for the typical time Tj in terms of a double integral over the 
frequen cy of the least-loaded c l ass which have been evaluated numeric ally for 



a = 1 (IStephan et all Il993l ; IGordo and Charlesworth! . l2000al ). Here 



we estimate these integrals analytically and find that for a broad range of pa- 
rameters, the average inter-click time is of a scaling form (129]) for any a > 1. 
For a = 1, it is shown that the ratchet speed is a function of the number 
n = Ne~ u l s in the least-loaded class and the selection coefficient s, and 
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not n alone as assumed in previous studies (IHaighi . Il978l ). With epistatic 
interactions, the time Tj is found to increase faster than any power law with 
J for any a > 1. Thus, an arbitrarily small amount of epistasis is sufficient 
to halt the ratchet with the ratchet speed approaching zero as ~ l/t for large 
time t. 



MODELS 



We consider a haploid asexual population evolving via mutation-selection 
dynamics starting with an initial condition in which all the individuals in 
the population have zero mutations. The genome length is assumed to be 
infinite so that back mutations can be ignored. If the population has a finite 
size N, it evolves stochastically following the discrete time Wright-Fisher 
dynamics. An offspring in generation t + 1 chooses a parent in the previous 
generation with a probability proportional to the fitness of the parent. Then 
the probability P(n) that a parent p carrying k mutations and with fitness 
W(k;p) has n descendants in one generation is given by 

(N\ { W(k;p) Y{, W(k;p) \ N - n ^ 
Pin)= [n)[-N{Wy) l^iVwi () 

where (W) = X^fcLo W(k)X(k, t) is the average fitness of the finite population 
in generation t. Here we have defined X(k,t) as the fraction of population 
with k mutations in a single sampling of Wright-Fisher process. From the 
above equation, it follows that the average number of offspring produced in 
one generation is proportional to the parent's fitness and the relative variance 
in offspring number decays as 1/N. This fact will be useful in defining the 
diffusion coefficient (1211) within diffusion approximation discussed in later 
section. Following replication, mutations are introduced where the number 
of new mutations is a random variable chosen from a Poisson distribution 
with mean U. In the simulations, the above process was implemented but the 
order of mutation and selection was reversed. An individual picked randomly 
from the population at time t was first mutated and the resulting mutant was 
allowed to survive at t + 1 with a probability equal to its fitness. This process 
was repeated until the generation t + 1 has N members and the population 
fraction X(k,t + 1) was recorded. It is useful to define Xj(k,t) = X(J + k,t) 
where J is the minimum number of mutations in the population at time t so 
that Xj(k, t) = for k < 0. If Xj(0, Tj) becomes zero, the least-loaded class 
J is lost and the ratchet has clicked at time Tj. 
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The ratchet effect due to which the least-loaded class is lost is essentially 
a stochastic problem arising due to the finite number N of individuals in 
the population. However as we shall describe later, the population fluctuates 
close to the deterministic frequency between two clicks of the ratchet. For 
this reason, we also study the problem of an infinite population for which the 
density fluctuations vanish and the average population fraction X(k ,t) with 
k mutations at time t obeys a deterministic quasispecies equation (jEiGEN . 
197ll : IJain and KkugL l2007h . Similar to the finite population problem, we 
define Xj(k,t) = X(J + k,t) where Xj(k,t) = for k < 0. Th en neglecting 
the back mutations for a genome of infinite length (lHlGGsl . ll994f ). the fraction 
Xj(k, t) evolves according to the following difference equation, 



Xj(j,t + 1) 



fc=0 



TJ3-k 

' U r—j^W(J + k)Xj(k,t) . 



(3) 



In this equation, the population fraction with k mutations replicates with fit- 
ness W(k) and accumulates further mutations which are Poisson distributed 
with a mean U . The average fitness Wj(t) = J2T=o W(J + k)Xj(k, t) in the 
denominator ensures that the number density is conserved. 



STEADY STATE OF THE QUASISPECIES MODEL 

In this section, we calculate the steady state population frequency Xj(k) 
in the error class J + k. Unlike for the multiplicative fitnes s case, the fre - 
quency Xj{k) depends on J for the epistatic fitness function (IHaighi . Il978l ). 
In particular, the fraction Xj(0) (later abbreviated as Xj) in the least-loaded 
class is expected to increase with J for a > 1 and decrease for a < 1. This 
can be explained by a simple argumen t which has also been used to under- 
stand the error threshold phenomenon (lElGENl . Il97ll ; IJain and KrugI . 120071 ) 
in which the fittest genomic sequence can get lost beyond a c ritical mutatio n 
rate in populations evolving on epistatic fitness landscapes (jWlEHH . Il997l ). 
Consider the ratio W( J + k)/W( J) ~ (1 - s ) akJ for J > k. For synergis- 
tic interactions, the error class J + k in the neighborhood of the least-loaded 
class has a fitness much worse than the fitness of class J rendering selec- 
tion effective in localising population in the class with J mutations. With 
increasing J, the selection pressure increases further. Thus we may expect 
the population frequency Xj(k) to peak around J and Xj to increase with 
J. On the other hand, in case of antagonistic epistasis (a < 1), the fitness 
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landscape is nearly neutral at large J so that the least-loade d sequence ca n 
be lost even in the deterministic limit (finite error threshold) (I WieheI . 119971 ) . 
In the steady state, the quasispecies equation (j3J) reduces to 



k=0 



TJ3-k 

' U J—^W{J + k)Xj{k) 



(4) 



where Wj is the average fitness in the steady state w hen the least-loaded class 



is J. The equation f or j — immediately shows that (IKimura and Maruyama 
19661 : IHaighL Il978h 

Wj = W(J)e~ u . (5) 



For j = 1 in 



we have 



UW{J)Xj 



W(J) -W(J + 1) 



Plugging this expression in the equation for j = 2, after some algebra we 
obtain 



*j(2) 



W(J) - W(J + 2 
Similarly, the fraction in the error class J + 3 is 



1 

2 + W(J) - W(J+1) 



U 3 W(J)Xj 



+ 



W{J + 1) 



W(J) -W(J + 3) L3! 2\(W(J)-W(J+1)) 



+ 



iy(J + 2) 



1 

777 + 



W(J + 1) 



W(J) - W(J + 2) V2! W(J) - W{J + 1 



From the expressions for Xj(k) for k = 2, 3 shown above, it is clear that 
in the weak selection limit s — > 0, the leading order contribution to Xj(k) 
comes from the last term. In general, we can write 



Xj(k) 



n 



W{J + m) 



U k W(J)Xj 

W(J + k) 11 W(J) - W(J + m) 

k k 



■ (j + my 



(6) 



7 



where the population Xj in the least-loaded class is determined using the 
normalisation condition Y2T=o'^j(^) = 1- 

Using the preceding equat ion for the multiplicative fitness function, we 
obtain the well known result (IKimura and MaruyamaI . Il966l ) 



-U/s 



k\ 



(7) 



The fraction Xj(k) for all k is seen to be independent of J (IHaighI . Il978l ). 
For synergistically epistatic fitness landscape with a = 2 in fitness function 
([[]), we have 



Xj(k) 



m=l 



u 



2Jm + m 2 \s J k\(2J + k)\ 



(2J)\Xj 



(8) 



On summing both sides over k, we find 



X] 




(9) 



wher e I n (z) is the modified Bessel function of the first kind (IAbramowitz and Stegun 



19641 ) . The fraction Xj(k) with J + k mutations is then given by 



■V, (k)= (- 
s 



J+k 



k\(2J+k): j 



2.1 



2J^ 



(10) 



and is plotted in Fig. [T^ as function of J + k at a fixed [// s. From the above 
equation, we find that for a given J, the fraction Xj(k) is centred about 
k\ = ^Jj 2 + Jl - J where J 2 = \/U/s. For J < J 2 , the distribution Xj(k) 
peaks about k* 2 w J 2 while for J ^> J 2 , it is maximised at J|/(2J) (see 
Fig. dh). Thus for large J, as argued at the beginning of this section, the 
distribution Xj(k) localises close to k = 0. The behavior of the least-loaded 
fraction Xj shown in Fig. [lb also depends on J 2 . For large J 2 (i.e. weak 
selection), the fraction Xj increases towards unity slower than for small J 2 . 
Using the asymptotic expansion of the Bessel function I n (z) for large orders 
(IAbramowitz and Steguni . ll964J ) in ([9]), we have 



X~ 



e 2Jy/lW 



2.1 



y 



(ii) 
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Figure 1: Steady state of the quasispecies model on synergistic fitness land- 
scape with a = 2: (a) Fraction Xj(k) as a function of J + k for U/s — 50 
given by ffTUj) (b) Fraction Afj of the least-loaded class J calculated using (jUJ) . 



where we have defined y = J 2 / 'J r . For J ^> J 2 , the above expression can 
be simplified to give Aj ~ exp(—U/(2sJ)) which asymptotically approaches 
unity. Thus with increasing J, most of the population tends to stay in the 
least-loaded class. 

For arbitrary a > 1, it does not seem possible to obtain explicit expression 
for Xj{k). However using the insights obtained from a = 2 case, we can find 
Xj for large J. We expect that for any a > 1, there exists a least-loaded 
class J a such that the population frequency Xj{k) with J ^> J a is nonzero 
for k <^ J. In such the denominator under the product sign in (J6]) 

can be expanded for m J to leading orders and yield 



As A'j(fc) decays fast with k, we can sum over both sides of the above solution 
to obtain 

Xj w exp [-[//(asJ^ 1 )] , J > J a . (13) 

This expression matches the exact results for a = 1 and 2 discussed above. 
The product in (jBJ) seems hard to calculate for J <^ J a . But for J = 0, we 
immediately have 

£A fc X 



Mk) « ( - ) ^ (14) 
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Figure 2: Time evolution of the fraction X(k, t) for k = 12 and 13 to illustrate 
the advance of the ratchet on epistatic fitness landscapes. Here a = 2, N = 
256, U = 0.25 and s = 0.005. 



which peaks at k* a equal to (U/s) l l a . By analogy with the a = 2 case, this 
suggests J a = (U/s) l l a . 

TIME BETWEEN SUCCESSIVE CLICKS OF THE RATCHET 



In this section, we first describe the process by which ratchet clicks and 
then calculate an expres sion for t he ty pical time between successive clicks 
using a diffusion theory (IEwensI . Il979l ). Let X(k,t) denote the fraction of 
population with k mutations at time t in a single realization of the Wright- 
Fisher process. Figure [2] shows the time evolution of X(k,t) for k = 12 and 
13 starting from an initial condition in which all the N individuals carry zero 
mutations. As Fig. [2^ illustrates, the fraction X(12,t) increases from zero 
to a steady state fraction about which it fluctuates until a time r n = 1194 
after which it relaxes to another steady state before finally dropping to zero 
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at Ti 2 = 2130 due to stochastic fluctuations. At t = r 12 , the ratchet is said 
to have clicked as the least-loaded class with 12 mutations gets irreversibly 
lost and the class with 13 mutations shown in Fig. [2b becomes the new least- 
loaded class which itself gets lost at T13 = 5235. Since the ratchet is clicking at 
a slow rate, X(k, t) has an opportunity to equilibrate. As Fig. [2b shows, soon 
after time r 12 , the fraction X(13,t) fluctuates about a steady state fraction 
which is close to the deterministic frequency Xi 3 given by Similarly, 
X(12, t) in Fig. [2^ oscillates about A12 after the 11th error class is lost until 
time T12. As a click of the ratchet is signaled by the change in the average 
value of X(k,t), it follows that the 11th class is lost at T\\. For m < t < Ti 2 
as there are at least 12 mutations in the population, the fraction X(13,t) 
equilibrates about the frequency Xi 2 (l). Thus the population fraction X(k, t) 
for fixed k passes through a series of steady states with frequency Xj(k — 
J), J < k before reaching the final absorbing state X(k,Tk) = 0. Note that 
this description of the mechanism by which the ratchet clicks assumes that 
the population NXj in the currently least-loaded class J far exceeds one and 
thus has a chance to attain equilibrium before the ne xt click. 



We will use the diffusion approximation proposed in lSTEPHAN et al\ (119931 ) 
to find the average inter-click time Tj = tj — tj_i between the ( J — l)th and 
Jth click of the ratchet where — stands for averaging over stochastic histo- 
ries. Let the random variable Xj G [0, 1] denote the population fraction in 
the least-loaded class J. If Xj = at time t, the current least-loaded class J 
is lost forever and the ratchet is said to have clicked at t. We are interested in 
calculating the average time T j required to reach the absorbing state Xj = 
starting from Xj at t' . The probability distri bution P ( Xj, t \X[ T , t') obeys the 



following backward Fokker-Planck equation (IRlSKEM Il996l ) 



dP(Xj,t\X f j,if) _ , ^dP(Xj,t\X'j,t') 



df iy J ' ' dx'j 

D2{X'jJ) d 2 P(Xj,t\X'j,t') 
+ 2 dXf (15) 



where 



D n (X'j,t') = Lt T ^ / dX. 

Jo 



(Xj - X'j) n P(Xj,1f + r\X'j,f) 



U T J X *' + ^- X ^W . (16) 
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As the coefficients D n in (IT51) are independent of t' (see below), the average 
inter-click time Tj defined as 



T(X>j 



dt' {-t!)P(0,0\X'j,l!) 



obeys the following ordinary differential equation, 



-1 = M 



dT D 2 {X'j) d 2 T 



dX'j 



+ 



dX'j 



(17) 



(18) 



Since X'j = is an absorbing state, the solution to the above equation is sub- 
jected to the boundary condition T(0) = 0. Furthermore as the population 
in the Jth class equilibrates about the mean Xj after the (J — l)th click, we 
can choose the initial distribution of random variable X'j to be S(X'j — Xj). 
Then the time T j during which J is the least loaded class obtained by solving 
(jTSjl is given by (IEwensI . Il979h 



Xj 



dY ij(Y) 



dZ 



ip(Z)D 2 (Z) 



(19) 



where ip(Y) = exp -2 J Y dX D 1 (X)/D 2 (X) 

We will now determine the coefficients Di and D 2 . The drift coefficient Di 
defined in ffTBT) measures the change in the average fraction of the least-loaded 
class over a generation. As the population is in local equilibrium, this can 
be determined using t he quasispecies equation ([3]) for j = 0. Thus the drift 



coeffic ient is given by (IStephan et all Il993l ; IGORDO and Charlesworthi . 
2000ah 



D l {X J = Xj(Q,t)) 



Xj(0,t+1)- 
~- u W(J) 



X 



J- 



Xj$,t) 
- Wj(t) 



(20) 



As expected, D\ vanishes when the population is either in the steady state 
(see ()5])) or in the absorbing state (Xj = 0). The equation (1201 for Di does 
not close in Xj but one can obtain an app roximate expression for D\ using 



the linear response theory (IRiskeni . Il996l ). As Di is proportional to the 



deviation from a steady state quantity, we can write (IStephan et all 119931 ) 



Dt oc e~ u W(J) - Wj(t) = C ( 1 



Xj_ 

Xj 



(21) 
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where C is a constant. Thus the drift coefficient can be written in terms of 
Xj as 

xy 



D^CXj^l-^-j . (22) 

As Fig.[2]shows, when the ratchet clicks (Xj = 0) and the Jth class is lost, the 
population quickly relaxes to the equilibrium frequency of the (J + l)th class 
so that the deviation in the fitness C ~ Wj-W J+1 = e~ u [W(J) - W(J + 1)] 
where we have used flS}. For s — > 0, expanding the fitness W(J) to leading 
orders in s, we get 

D 1 {X J )^AXj[l-^ (23) 

where A — cs [(J + l) a — J a ] ~ acs J"" 1 for large J. 

The diffusion coefficient D 2 in (fTBT) gives the fluctuations in the frequency 
of the least-loaded class about the mean value. These fluctuations arising due 
to the finiteness of the population can be determined using (j2J) which gives the 



varia nce in the number of offspring produced in one generation as (jEwENS . 



19791 ) 

D 2 (Xj) = Xjil N Xj) « § . . (24) 

The last expression on the right hand side of the above equation captures 
the fact that the fluctuations vanish when either the population size N is 
infinite or the population is in the absorbing state Xj = 0. 

Using the coefficients (1231) and (|24l) in (fT9l) . the average inter-click time 
Tj can be written as 

' Xj „, f l dZ 



dY ^(Y) ^ a — r\Z) (25) 



where 



and 



ip(Y) = exp 



p\--- 

A. j A. j 



(26) 



13 = NXjA = NXjcs [{J+ l) a - J a ] . (27) 



In the absence of epistasis, both A and Xj are independent of J so that 
typical time spent between any two successive clicks is constant and the 
ratchet turns with a finite speed equal to 1/T. For epistatic fitness a / 1, 
Tj depends explicitly on J and is expected to increase with J for a > 1 
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while decrease for a < 1 (IKondrashovI . Il994l ). In the following discussion, 
we will restrict ourselves to a > 1. 

After some simple manipulations, we can rewrite fl25|) as 



2NXj / dYe^ 2 



-i 



Y 



dZ 
l + Z 



-pz 2 



(28) 



which implies that the scaled time Tj/(NXj) is a function of two variables 
namely (3 = NXjA and 5 = (1 — Xj)/Xj. The nature of Tj depends on the 
parameter (35 2 which can be seen as follows. Consider the Gaussian e _/3Z in 
the rightmost integral in (128]) which is centred about Z = and has a width 
If the upper limit 5 of this integral exceeds the width i.e. [35 2 3> 1, 
the integral can be cutoff at l/Vft thus eliminating the dependence on 5. In 
such a case, Tj is of the following scaling form, 



Tj^NXj F{NXjA) 



(29) 



where F(f3) is the scaling function determined below. If (35 2 <C 1, then the 
inter-click time depends on 5 as well. Since f35 2 ~ N sJ ' 1 / Xj and Xj is 
bounded above by one, j35 2 3> 1 for large J when epistatic interactions are 
synergistic. For a = 1, the parameter j38 2 exceeds unity if Ns ^> 1. Here we 
will restrict ourselves to the /35 2 ^> 1 case which has nice scaling properties 



although the double integral in (|28|) can be estimated for f35 2 -C 1 also. 

We now proceed to find the scaling function J-'(P). If the width of the 
Gaussian is large i.e. (3 1, we can approximate e _/3Z ~ 1 for Z <^ l/y/j3 
in the rightmost integral in ( |28l) . Since 5 ^> l/y/j3, as argued above, this 
integral needs to be carried out from Y to 1/ \fj3. This yields 



Tj 



2NX T / dYe^ 2 



In 1 + 



1 

7^ 



-in(i-y) 



2NX, 



NXj(2 



1 + In 1 + 



1 

ln/3) , [3 < 1 



(30) 



To find the scaling function in the opposite limit (3 3> 1, we first consider the 
inner integral in (|28|) . 



dZ 
l + Z 



-0Z* 




Z 

VP. 



i 



(3) - erf (YW/3 



(31) 
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where erf(z) is the error function (IAbramowitz and Steguni . Il964j ) and 
we have kept terms to leading orders in /3<5 2 (3> 1) and Y for @ ^> 1. On 
using the above integral in fl28l) . we obtain 



T., 



NXj^e 13 f p dYe~ Y 



erf <V/3 +erf (y) 




erf + erf ( y/0 - Y 



(32) 



where we have used that both (38 2 and j3 are large. From the above discussion, 
we find that the scaling function 



HP) 



2 - In/? , /?<1 



(33) 



is a U-shaped function of /3 reaching a minimum when /3 ~ 1. We now 
discuss the above results in more detail for a = 1 and a > 1. 
Non-epistatic fitness landscapes: As we have already argued, when 
a = 1 the ratchet clicks with a finite speed equal to 1/T. When the rate 
at which the ratchet clicks is la rge, analytical r e sults can be obtained us- 
ing the traveling wave approach (IRouziNE et all 120031 ) . For slowly clicking 
ratchet which is the subject of this article, the prob lem was formula t ed an - 
alytically within a diffusion approximation first by IStephan et all ( 119931 ) . 
However a better agree ment between the diffusion theory and t he simulation 
results was obtained in IGordo and CharlesworthI ( 2000al lbl). A possi- 
ble reason for this difference is that IStephan et all ( 119931 ) included terms 
besides those in ([22]) in the expansion of the drift coefficient which is not 
consis tent with the assumption of linear resp onse, while the expression for 
D 1 in IGordo and CharlesworthI (j2000al ) is the same as ([22]). In fact, 
the expression for t he inter-click time given by (|28l) with c = 0.6 is identical 
to that reported in lGoRDO and CharlesworthI ( ]2000al ). However the in- 
tegrals were computed numerically by these authors while here we estimate 
them analytically and find that the average inter-click time is given by 



n (2-ln/3) , j3 < 1 
v^Fno /r 3 /V , (3 > 1 



(34) 
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Figure 3: Data collapse of the scaled time T/n when plotted against the 
scaling parameter (3 = hqcs for the multiplicative fitness landscape where 
c = 0.6. The simulation data has been averaged over 100 stochastic runs. 



where no = NXj = Ne~ u ' s is the number of individuals in the least-loaded 
class and (3 = uqcs is the scaling parameter. The inter-click time calculated 
numerically using the integral ff28l) and the above expression flM|) is shown 
in Tableland the two are seen to be in good agreement. 

Our first result concerning the Muller's ratchet is the scaling form for 
time T when parameters N, U and s are chosen such that (35 2 ~ Ns 3> 1 is 
satisfied. The results of our numerical simulations testing this scaling form 
are shown in Fig. [3] where we have differentiated data points at fixed (3 for 
clarity. The scaled time T/n indee d shows a very good dat a colla pse and 
a non-monotonic dependence on (3. IHiggs and Woodcock! (119951 ) unsuc- 
cessfully attempted to obtain a scaling form for the ratchet rate by using Ns 
as the scaling parameter for Ns ^> 1. The scaling function in (I34p however 
shows that the scaling parameter is a function of the population number uq 
in the least-loaded cl ass, and not the total population N. 



In many studies (IHaighI . Il978l ; IBellI . Il988l ; IGesslerI . 119951 ). the size 



n of the least-loaded class has been regarded as an important parameter in 
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Table 1: Comparison of simulation data plotted in Fig. [3] with analytical 
results. Here T sim refers to the inter-click time obtained in simulations, T int 
by numerically evaluating the integrals in f[2"51) for a — 1 and T sca using the 
scaling form (1341) . 



determining the ratchet speed. If n ^> 1, the population is close to the de- 
terministic limit and the ratchet clicks slowly whereas for no <C 1, the ratchet 
speed is high. However the simulations show that t he size nn of the least- 
loaded class is not sufficient to predict the ratchet rate (IStephan et a/.LIl993 ; 



Gordo and Charlesworth 



2000al ). This is indeed captured by diffusion 



approximation as (1341) is not a function of no alone. However, if s is kept fixed, 
T increases monotoni cally with np in accordance with the above expecta- 
tion a nd simulations ([Gabriel et all Il993l ; IGordo and Charlesworth! . 



2000a|) . The dependence of T on s for given nn is however non-monotonic sim - 



ilar to that seen in numerical studies (IGordo and Charlesworth . l2000al ) . 
To understand this behavior qualitatively, consider the situation where the 
population number uq is kept constant by keeping N and U / s fixed. As in- 
creasing s tends to localise population and hence increase T, increasing U has 
the opposite delocalising effect which decreases T. At a given U and s, one of 
these two competing forces wins. According to ( l34l) . as the scaling function 
overturns when nos ~ 1, the mutation takes over for U > U* = sln(iVs) 
whereas below U*, selection dominates and T is large. 
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The solution (1341) gives an initial logarithmically slow drop in s and an 
exponential increase for larger s with the minimum of the U-shaped curve 
oc curing at a selection coefficient which sc ales as I/tlq. The simulation results 
of IGordo and Charlesworth feOQOal ) (also see Tabled]) however show a 
much faster drop at small s. A good agre e ment w ith simulation data was 
obtained in IGordo and Charles worth ( l2000bl ) by adding T and time 
T a ~ 1/s required by the population to relax to new steady state just after a 
click. However, the full expression for T a is not of scaling form (|29|) although 
the simulation data in Fig. [3] shows an excellent data collapse even for small 
s. In view of this, a better understanding of the time T a is desirable. Of 
course, for both n and s fixed, the time T is predicted to be independent of 
N and U which is confirmed by simulations (IGordo and Charlesworthi . 



2000a|). 



Epistatic fitness landscape s: For synergistic interactions, the ratchet is 
expec ted to halt at large times (|Charlesworth et a/lll993l ; lKoNDRASHOvl . 



1994J ) . FigureHk shows the population fraction X(k, t) with k deleterious mu- 



tations at several time slices. Two points are noteworthy: the ratchet does 
not turn with a constant speed as is evident by the rate at which the pop- 
ulation accumulates average number of mutations at late times . Secondly, 
unlike for the multiplicative fitness case (IRouziNE et all 120031 ). the popu- 
lation fraction does not maintain its shape as the width of the distribution 
X(k,t) decreases with increasing time. Thus the number frequency of an 
asexual population under epistatic selection does not behave like a traveling 
wave moving with a constant speed. 

The simulation data for average inter-click time shown in Fig. Hp increases 
with the minimum number J of deleterious mutations in the population thus 
indicating the arrest of the ratchet at large times. The time Tj obtained by 
integrating ( 1281) for a = 2 is also shown for comparison and we find that it 
agrees well with the simulation results. As the scaling parameter f3 defined 
in (1271) increases as a power law with J for J ^> 1 , we have (3 3> 1 for large J. 
In such Tj ~ NXj(3-*/ 2 eP which increases exponentially fast with (3. 

As Xj — > 1 for large J, the inter-click time T j ~ e Nacs,j a increases faster 
than any power law with J and for any a > 1. Thus an arbitrarily small 
a — 1 is capable of slowing down the ratchet under synergistic selection. 

To estimate how the ratchet speed approaches zero, we use the following 
argument. The average speed v = d/t where d is the number of minimum 
mutations accumulated in time interval t. Since Tj is the time between the 
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Figure 4: Loss of least-loaded class on epistatic fitness landscapes with a = 
2, iV = 256, U = 0.25, s = 0.005: (a) Average fraction X(k, t) in the jfeth class 
at t = 50, 100, 200, 1250, 2500, 5000, 10000 (from left to right). The data have 
been averaged over 1000 stochastic histories, (b) Average inter-click time Tj 
as a function of J where each data point has been averaged over 500 histories. 
Comparison with Tj obtained by evaluating integral (1251) numerically and 
using the scaling form (]29|) at large J is also shown. A reasonably good 
agreement is obtained with c = 0.9. 
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(J — l)th and Jth click of the ratchet, we have 



d 

where d is fixed. Assuming the distribution for click times has nice scaling 
properties, we may write 

" (36) 



Tn + ...+T ( 



If Tj does not depend on the least-loaded class J, we have v — 1/T as 
expected for a = 1 case. For a > 1 as the average inter-click time increases 
faster than power law with J, the sum T + ... + T d ~ d 2 - a e Nacsda which 
is reasonable as the sum is dominated by Td for large d. This gives v ~ 

gp-l e -Nacsd a Qr - n ^ erm g Q f + 



- 1 / 

f ~ - - 



lnt \ 



l/(a-l) 



£ \iVacs 



J • (37) 



Thus in the presence of epistasis, the average speed of the ratchet approaches 
zero as 1/t with ^-dependent logarithmic corrections. 

DISCUSSION 

In this article, we considered the effect of drift and epistasis on the loss 
of the least-loaded (or the fittest) class in an asexual population. When 
the population size is infinite, the drift is absent and the population evolves 
due to the elementary processes of selection and mutation. As selection 
tends to localise the population at the fittest sequence while mutation has 
the opposite tendency to delocalise it, an error threshold ma y exist beyon d 



which the fittest class cannot be sustained in the population (lElGENl . Il97ll ). 



Such a phase transition is known to occur for asexual populations evolving 



deterministically on fitness landscapes defined by (CQ) when a < 1 (IWiehe . 



19971 ). However for a > 1, the population frequency Xj in the least-loaded 
class J remains nonzero for any finite U, s. In fact, for synergistic interactions, 
the frequency Xj given by ffT3l) increases with J towards unity. 

If however the population size is finite, as illustrated in Fig. [2J the pop- 
ulation frequency Xj(t) fluctuates with time and can become zero even for 
a > 1. Numerical simulations have shown that this loss occurs at a constant 
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speed when a = 1 (IHaighi. Il978l) but at a deceler ating rate when a > 1 



(ICharlesworth et all Il993l ; IKondrashovI . Il994j ) . In this paper, we fo- 



cused on the stochastic dynamics of the loss of the fittest class for a > 1 
and calculated explicit analytical expressions for the typical time Tj dur- 
ing w hich the least-loaded class J survives using a diffusion theory (IEwens . 



19791 ). Although this appr oach has been consider ed previously to attack the 



Muller's ratchet problem (IStephan et all Il993l ). the resulting solution in 



the form of a double integral was evaluated numerically which does not allow 
one to infer the functional dependence of Tj on parameters N, U, s and J. 

When the interactions between gene loci are assumed to be absent and 
Ns ^> 1, the inter-click time T is found to be of the scaling form (1M|) . 
Although this result is derived using diffusion theory which is based on sev- 
eral approximations, the numerical simulations show an excellent data col- 
lapse suggesting that ff29l may be an exact statement. For fixed n , the 
time T is seen to be a U-shaped function of s arising due to competition 
between mutation and selection. Such a behavior is reminiscent of error 
threshold phenomenon in infinite populations discussed above. Although 
the least-loaded cl ass is never lost in the determ inistic limit on multiplicative 
fitness landscapes ( IWagner and KrallI . 119931 ). the selection- mutation com- 
petition manifests itself in the time duration during which the least-loaded 
(fittest) class can support a finite population. For given s, the survival time 
T initially increases linearly with uq approaching the deterministic limit of 
N — > oo with an exponential rise with uq as increasing N decreases the effect 
of drift. 

Muller's ratchet has been proposed as a possible mechanism for the de- 
generation and eventual extinction of asexual organisms and nonrecombining 
parts of sexually rep roducing popu l ation s. For bacterial populations with 
N ~ 10 6 , U » 0.003 (IDrake et all Il998h and s « 10~ 3 , using (JMD we find 
it takes 10 15 generations for one click to occur which does not seem plausible. 
However, reducing s by a factor half gives T ~ 5000 and further reduction 
to s ~ 0.25 x 10 -3 gives just 50 generations for a single click. Muller's 
ratchet has also been invoked to understand the degeneration of the nonre- 
combining neo-Y chr omosome which origina t ed about a million years ago in 
Drosovhila Miranda (|Charlesworthi . I1978| ; IGordo and Charlesworthi . 
2000al ). Assuming that a few thousand deleterious mutations have occurred 



over this time span, the time between successive turns of the ratchet is of 
the order of 10 3 generations. From (|34p . the time T ~ e n ° cs for large T 
which gives uqs ~ 14. If s ~ 10~ 2 , the population size iV required for the 
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Muller's ratchet to be a viable mechanism works out to be ~ 700e 100U which 
depends sensitively on U due to the exponential dependence. For instance, 
for U = 0.07, the required population is of the order 5 x 10 5 while it reduces 
by a factor twenty for U = 0.04. A similar sensitive dependence of extinction 
time on s for given N and U has been no ted in the prob lem of the degenera- 
tion of human mitochondrial DNA also (ILoeweI . 120061 ) . This suggests that 
very precise estimates of U and s may be required to determine whether the 
Muller's ratchet might be in operation. 

The scaling form (|29j) holds for a > 1 also but the speed of the ratchet 
under epistatic selection is found to decay rapidly with time for any a > 1. 
Although the Muller's ratchet under synergistically epistatic selection has the 
interesting feature of arresting the loss of least-loaded class, the generality 
of this mechanism seems unclear as gene ral support for synergistic epi stasis 
has not been found in the experiments (IDE VlSSER and ElenaI . 120071 ) and 
the slowing down effect is sensitive to the inclusi on of biologically relevant 
details such as distribution of mutational effects (IButcherI . Il995l ). Recent 
experimental evidence suggests that fitne ss decline down t o a p lateau can 



be attributed to the presence of epistasis (ISilander et all 120071 ) . This can 



be due to negative epistasi s (jKoNDRASHOvt I1994T ) or compensatory epista- 
sis (ISilander et all 120071 ) . It would be interesting to compare the results 
discussed in this work with models that inc l ude c ompensatory mutations 



(IWagner and Gabriel! . Il990l ; IWilke et all 120031 ). 
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